## Helper function to shorten amount of code required to specify MxModels in OpenMx
JmxModel <- function(..., manifestVars, latentVars, data, dtype = c("raw", "cov", "cor"),
  type = "RAM", model = NA, name = NA, run = TRUE, intervals = FALSE) {
  tmp.args <- list(...)
  new.args <- unlist(tmp.args, recursive = FALSE)

  ## handle manifest variables
  m.index <- grep("manifestVars", names(new.args))
  if (length(m.index) > 0) {
    MVs <- c(unlist(new.args[m.index], recursive = TRUE))
    attributes(MVs) <- NULL
    new.args[m.index] <- NULL
    new.args$manifestVars <- MVs
  }

  if (!missing(manifestVars)) {
    if (any(grepl("manifestVars", names(new.args)))) {
      manifestVars <- c(new.args$manifestVars, manifestVars)
    }
    new.args$manifestVars <- manifestVars
  }

  ## handle latent variables
  l.index <- grep("latentVars", names(new.args))
  if (length(l.index) > 0) {
    LVs <- c(unlist(new.args[l.index], recursive = TRUE))
    attributes(LVs) <- NULL
    new.args[l.index] <- NULL
    new.args$latentVars <- LVs
  }

  if (!missing(latentVars)) {
    if (any(grepl("latentVars", names(new.args)))) {
      latentVars <- c(new.args$latentVars, latentVars)
    }
    new.args$latentVars <- latentVars
  }

  ## Default to RAM model if no type specified
  if (!"type" %in% names(new.args)) {
    new.args$type <- type
  }
  if (all(!c("model", "name") %in% names(new.args))) {
    new.args$model <- model
    new.args$name <- name
  }

  ## if the data type (dtype) is not specified AND there are no mean paths
  ## use the raw data and set default paths from one to all manifest variables
  if (missing(dtype) && !any(sapply(new.args, function(i) {
    if (inherits(i, "MxPath")) "one" %in% i@from else FALSE}))) {
    new.args <- c(new.args, mxPath(from = "one", to = new.args$manifestVars,
      labels = paste("m", new.args$manifestVars, sep = '')))
  }

  ## if no MxData object is supplied create one
  ## defaulting to the raw manifestVars in data
  if (!any(sapply(new.args, function(i) inherits(i, "MxData")))) {
    dtype <- match.arg(dtype)
    data <- data[, new.args$manifestVars]
    if (dtype %in% c("cov", "cor") && any(is.na(data))) {
      warning("Listwise deletion is being used for missing data. Use type = 'raw' to use FIML")
    }
    new.args <- c(new.args, switch(dtype,
      raw = mxData(observed = data, type = "raw"),
      cov = mxData(observed = cov(na.omit(data)), type = "cov", numObs = nrow(na.omit(data))),
      cor = mxData(observed = cor(na.omit(data)), type = "cor", numObs = nrow(na.omit(data)))
    ))
  }

  M <- do.call("mxModel", new.args)
  if (run) M <- mxRun(M, intervals = intervals)
  return(M)
}

saturateMe <- function(object) {
  mvars <- object@manifestVars
  object@data[, mvars]
  mn <- object@name
  satMod <- mxModel(paste("Saturated", mn), type = "RAM", manifestVars = mvars,
    mxPath(from = mvars, arrows = 2, lbound = 0L),
    mxPath(from = combn(mvars, 2)[1, ], to = combn(mvars, 2)[2, ]),
    mxPath(from = "one", to = mvars))
  satModFit <- mxRun(satMod)
  return(satModFit)
}


## Define a function to calculate the wald statistic
## optionally creates boxplots and a QQ plot to examine normallity

WaldTest <- function(y, x, h, a = 0, data, plot = FALSE, showparams = FALSE) {

  if (is.language(y)) {
    model <- terms.formula(y)
    vars <- attr(model, "variables")
    vnames <- sapply(vars, deparse)[-1L]

    if (missing(data)) {dat <- eval(vars, .GlobalEnv)
    } else {dat <- eval(vars, data)}

    dat <- as.data.frame(dat)
    colnames(dat) <- vnames
    y <- dat[, attr(model, "response")]
    cells <- interaction(rev(dat[, -attr(model, "response")]))
    x <- sapply(levels(cells), function(x) ifelse(cells == x, 1, 0))
  }

  e <- environment()
  test <- sapply(c("y", "x"), function(x) {
    (is.vector(get(x, e)) || is.matrix(get(x, e))) &
    is.numeric(get(x, e))})
  if (any(!test)) {
    stop("Please check that '",
         paste(names(test)[!test], collapse = "' & '"),
         "' is/are numeric vector(s)/matrice(s).")
  }

  n <- length(y); k <- ncol(x)
  gnames <- paste("Group", 1:k, sep = '')

  xx <- ginv(crossprod(x))
  mu <- xx %*% t(x) %*% y
  dimnames(mu) <- list(gnames, "Mean")
  dev <- y - x %*% mu
  vareps <- as.vector(crossprod(dev)/(n - k))

  WaldFit <- function(h, a, mu, xx, df2 = n - k, var) {
    df1 <- nrow(h)
    hfit <- h %*% mu - a
    statistic <- (t(hfit) %*% ginv(h %*% xx %*% t(h))
      %*% hfit) / (df1 * vareps)
    pval <- pf(statistic, df1 = df1, df2 = df2,
      lower.tail = FALSE)

    output <- matrix(c(statistic, df1, df2, pval), nrow = 1,
      dimnames = list(NULL, c("F", "df1", "df2", "p.value")))

    return(list(output = output, h = h, a = a))
  }

  if (is.list(h) && is.list(a) && identical(length(h), length(a))) {
    tmp <- vector("list", length(h))
    for(i in seq_along(h)) {
      tmp[[i]] <- WaldFit(h[[i]], a[[i]], mu, xx, n - k, vareps)
    }
  } else {tmp <- list(WaldFit(h, a, mu, xx, n - k, vareps))}

  hlabels <- paste("H", seq_along(tmp), sep = "")
  table.out <- do.call("rbind", lapply(tmp, `[[`, "output"))
  rownames(table.out) <- hlabels
  h.out <- lapply(tmp, `[[`, "h")
  names(h.out) <- hlabels
  a.out <- lapply(tmp, `[[`, "a")
  names(a.out) <- hlabels

  output <- list(summary = table.out, h = h.out, a = a.out,
    mu = mu, variance = vareps, y = y, x = x)

  if (plot) {
    old.par <- par(no.readonly = TRUE)
    on.exit(par(old.par))
    par("ask" = TRUE)
    group <- factor(as.vector(x %*% 1:k),
      levels = 1:k, labels = gnames)
    plot(x = group, y = y, main = "Boxplot(s) of Scores")
    qqnorm(y = dev)
  }

  cat("Summary Table \n\n")
  print(round(output$summary, 5))
  if (showparams) {
    cat("\n Variance: ", round(output$variance, 5))
    cat("\nMu matrix: \n")
    print(output$mu)
  }

  return(invisible(output))
}

## Function to compare timings of different functions on the same data
compare <- function(data, functions) {
  myenv <- environment()
  x <- length(functions)
  out <- vector(mode = "list", length = x)
  times <- vector(mode = "list", length = x)
  names(times) <- functions
  for(i in 1:x) {
    FUN <- match.fun(functions[i])
    times[[i]] <- system.time(out[[i]] <- FUN(data))
  }
  checks <- expand.grid(1:x, 1:x)
  cat(all.equal(out[checks[, 1]], out[checks[, 2]]), fill = TRUE)
  return(times)
}

## Function to aid in converting data into a TeX table
TeXPrint <- function(dat, digits = FALSE) {
  f <- function(x) {
    paste(x, collapse = " & ")
  }

  if (digits) {
    dat <- as.data.frame(dat, stringsAsFactors = FALSE)
    dat <- as.data.frame(lapply(X = dat, FUN = function(x) {
      if (isTRUE(is.numeric(x))) {
        round(x, digits = digits)
      } else x
    }))
  }

  tmp <- apply(X = dat, MARGIN = 1, FUN = f)
  cat(tmp, sep = " \\\\ \n")
}

## T contrast function##
## raw data works properly now, I think

ContrastT <- function(x=NA, lambda, m=NULL, s=NULL, n=NULL, raw=TRUE) {
  if(identical(raw, TRUE)) {
    xok <- list(NULL)
    for(i in seq_along(x)) {xok[[i]] <- x[[i]][!is.na(x[[i]])]}
    for(i in seq_along(xok)) {m[i] <- mean(xok[[i]])}
    for(i in seq_along(xok)) {s[i] <- sd(xok[[i]])}
    for(i in seq_along(xok)) {n[i] <- length(xok[[i]])}}
  df <- sum(n - 1)
  effect <- sum(m * lambda)
  s2.pooled <- weighted.mean(x = s^2, w = n - 1)
  sample.correction <- sum((lambda^2) / n)
  variability <- sqrt(sample.correction * s2.pooled)
  t.score <- effect / variability
  p.score <- pt(q = t.score, df = df, lower.tail = FALSE)
  r.score <- t.score / sqrt((t.score^2) + df)
  value <- list(t.score, p.score, r.score, s2.pooled, df)
  names(value) <- c("t.contrast", "p.value", "r.contrast", "pooled.variance", "df")
  return(value)
}

## Mean Polishing Function

## x is a matrix (or 2-way array).  If inverse=FALSE, the default,
## x should be the raw data and will be mean polished.  If inverse=TRUE,
## x should be the residuals, plus the row and column effects at the
## edges of the table.  If it is going to 'unmean' polish, it expects
## the same input as it would output for regular mean polishing.
## Accordingly the following should return the original matrix:
## mean.polish(x=mean.polish(x=matrix(1:12, ncol=3), inverse=TRUE)

MeanPolish <- function(x, inverse = FALSE, na.rm = TRUE) {
  if(isTRUE(inverse)) {
    gmean <- x[nrow(x), ncol(x)]
    rmean <- x[-nrow(x) , ncol(x)]
    cmean <- x[nrow(x), -ncol(x)]
    fun <- "+"
    dat <- x[-nrow(x), -ncol(x)] + gmean
  } else {
    gmean <- mean(x)
    rmean <- rowMeans(x)-gmean
    cmean <- colMeans(x)-gmean
    fun <- "-"
    dat <- x - gmean
  }
  dat <- sweep(x = dat, MARGIN = 1, STAT = rmean, FUN = fun)
  dat <- sweep(x = dat, MARGIN = 2, STAT = cmean, FUN = fun)
  if(isTRUE(inverse)) {return(dat)}
  dat <- cbind(dat, "RowEffects" = rmean)
  dat <- rbind(dat, "ColEffects" = c(cmean, gmean))
  return(dat)
}

########## Unexported functions ##########

## T test function ##

TStat <- function(m, s, n, two.sample = FALSE, digits = 4) {
  stopifnot(is.numeric(m) && is.numeric(s) && is.numeric(n),
            length(m) == length(s) && length(m) == length(n),
            (isTRUE(two.sample) && length(m) == 2) ||
              identical(two.sample, FALSE))
  if (identical(two.sample, FALSE)) {
    variability <- sqrt((s^2)/n)
    df <- n - 1
    t.score <- m / variability
    ltail <- ifelse(t.score >= 0, FALSE, TRUE)
    p.score <- pt(q = t.score, df = df, lower.tail = ltail)
    t.crit <- qt(p =.05/2, df = df, lower.tail = FALSE)
    error <- t.crit * variability
    LL <- m - error
    UL <- m + error
  } else {
    variability <- sqrt(sum(weighted.mean(x = s^2, w = n - 1) / n))
    df <- sum(n) - 2
    t.score <- (m[1] - m[2]) / variability
    ltail <- ifelse(t.score >= 0, FALSE, TRUE)
    p.score <- pt(q = t.score, df = df, lower.tail = ltail)
    t.crit <- qt(p = .05/2, df = df, lower.tail = FALSE)
    error <- t.crit * variability
    LL <- m[1] - m[2] - error
    UL <- m[1] - m[2] + error
  }
  r.score <- t.score / sqrt((t.score^2) + df)
  CI <- paste("[", round(LL, digits = digits), ", ",
              round(UL, digits = digits), "]", sep = '')
  value <- list(t.score, p.score, CI, r.score)
  names(value) <- c("t.value", "p.value", "c.interval", "r.effectsize")
  return(value)
}

## r Confidence Interval Function ##
RCI <- function(r, N, conf.level = .95, digits = 4) {
  zr <- atanh(r)
  error <- qnorm((1 - conf.level)/2, lower.tail = FALSE) / sqrt(N - 3)
  CI <- tanh(zr + c(-error, error))
  cat("[", paste(round(CI, digits = digits), collapse = ", "),
    "]", sep = '', fill = TRUE)
  return(invisible(CI))
}

##Function that can either load or install then load a package##

RequireInstall <- function(x) {
##  old.repos <- getOption("repos")
##  on.exit(options(repos = old.repos))
##  new.repos <- old.repos
##  new.repos["CRAN"] <- "http://cran.stat.ucla.edu"
##  options(repos = new.repos)
  x <- as.character(substitute(x))
  if(isTRUE(x %in% .packages(all.available=TRUE))) {
    eval(parse(text=paste("require(", x, ")", sep="")))
  } else {
    ## recommended before installing so that dependencies are updated
    update.packages()
    eval(parse(text = paste("install.packages('", x, "')", sep='')))
    eval(parse(text = paste("require(", x, ")", sep='')))
  }
}

## Function to simplify converting p-values to asterisks
star <- function(x) {
  symnum(x, legend = FALSE, na = FALSE,
    cutpoints = c(0, 0.001, 0.01, 0.05, 1),
    symbols = c("***", "**", "*", ""))
}

## Calculate F and p-value from the R^2
f.r2 <- function(r2, numdf, dendf) {
  F <- (dendf/numdf) * (-r2/(r2 - 1))
  p <- pf(F, df1 = numdf, df2 = dendf, lower.tail = FALSE)
  c(F = F[[1]], numdf, dendf, p = p[[1]])
}
